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ABSTRACT 

A  bubble  transpart  equation  vjas  deve'loped'  and-a^  formal' solution 
obtained  for  tbe  bubble  distribution-,  function. .  The  distribution- 
function  obtained  depends-  upon  the-  model  used_  for:  bubble  drag ^  -gas^ 
diffusion,  and  ocean  circulation . .  As  a-- specific  application-,  .  solutions 
were  obtained  for  a  one— dimensional,,  quiescent,,  steady  state  ocean 
with  a  bottom  source  a-t  197"  meters..  The:  char.acte-ristic:  equations ' 
were  integrated  n.umericaliy  to.  tssir  the.  accura^ey  of:'  approximate  - 
analytical  solutions..  The  calcuiated  bubble  densities  were  compared  : 
with  experimental  data  available,  for  this  problem..  When  gas  diffusion 
is  neglected,  the  "t±.earetical  and  experimen-tai  distributions  are.  in- 
general  agreement.. 
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I.   IITTEQDUCTTON 

Studies  have  been  made  by  many  res.earchei:?s  regHrdihg^  the  mechanisms 
of  gas  bubble  dynamics  in  a  fluid  from  both,  the:  the:or^1rical.  [Ref . .  1]^ 
and  the  experimental  [Refs.  2,   3^  and  M-I  aspect^..  The' meehanism^s  of ' 
gas  bubble  formation,  transport,,  and  extinction  are  iinportant:  as  gas 
bubbles  can  effect  the  propagation  of  sound  in.  the  ocean',  by  the  under-- 
water  scattering  of  sound  beams..  TheorelricaX  studies"  have  explained 
how  a  bubble  varies  in  size  as  it  rises  in.  a  quiescent  fluid,.  The 
phenomena  involved  includes  shrinkage  from,  gas  diffusion,  versus 
expansion  from  decreets ing  hydrostatic  pressure  and  shrinkage  from 
surface  tension  for  very  small  bubbles..  However,,  there  does  not. appear 
to  be  at  present  any  theory  explaining  hcn-7' the  distributiDn. of:  bubbles 
generated  by  a  bottom  or  distributed  source  in. the  ocean. can  be  expected 
to  vary  vrith  depth  or  initial  bubble  radius.. 

This  paper  outlines  the  theory  far   determining  theoretically  the 
distributions  of  bubbles  generated  by  a  bottom  source  in  the  ocean 
[Ref.  5].   First,  a  linear  biibble  transport  equation  is  developed  and 
then  it  is  solved  by  the  method  of  characterislrics...  A  model  for.  a  one-- 
dimensional  steady  state  ocean  is  specified  for  application  of  this 
theory.   Depending  on  the  assumptions  made  regarding  the  dynamics 
involved  in  the  bubble  motion^  such  as  the  importance  of  gas  diffusion 
or  surface  tension,  various  semi— analytical  solutions,  can  be  obtained 
for  the  distribution  of  bubbles  from  a  bottom  s^ource..  These.^  distri-- 
butions  are  then  compared  to  available  experimental  data"..  Ih"  particular, 
the  distributions  from  a  bottom  source  located  at  197  meters  were  calcu- 
lated using  various  approximations  and  the  results  were  compared  with 
the  experimental  data  in  Ref..  4.. 
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II.   THE  BUBBLE  TPAirSPORT  EQUATTON 

The  transport  equation  is  basicaULy  a  balance  equation'  [Refs..  5 
and  7].   It  was  desired  to  develop  a  transport  equation- for,  the. 
description  of  the  distribution  of  bubbles  as  a  funclion  of/ position, . 
velocity,  size  (bubble  radius)  and  time  in  a  circulalring"  field  that 
may  contain  sources  or  sinks.   A  distribution  function  that  describes 
the  bubble  population  can  be  defin^ed  and  If  the  bubbles^  satisfy  the 
following  general  assumptions,  this  distribution-  wil  1  satisfy  as. 
Boltzmann-type  transport  eq^uation  r 

(1)  The  bubbles  are  completely  independent  of  one:  another/, . 

(2)  If  bubble  interaction  occurs^  no  more  than- two  bubbles  will 
interact  at  any  one  time. 

Since  bubbles  in  the  ocean  exact  no  force  on  one  another-  and-,  the 
history  of  one  bubble  does  not  affect  the  history  of  another  bubble, 
the  two  above  assumptions  hold  and  a  transport  equation  was  found. . 

If  i(>(r,v,£  ,t  )d3rd3vd£  is  the  number  of  bubbles  at  time  t- with 
radius  d£  about  £,  velocity  d^v  about  v,.  and  position  d-^r  about  r, 
then  at  time  t'  e  t  +  dt,  ip(r'' ^v*,]?,'' ^t ' )  d3r''d3v'd£'  is  the  number-  of 
bubbles  at  time  t  +  dt  with  radius  d£ '  about  £'  ,.  velocity  d^v'  about  v' 
and  position  d^r*  about  r'',  where  r'^  v',.  and  £*'  are:: 

v^    =  V  +   vdt 

v'  =  V  +.  adt  (1) 

I*    =  g  +   vdt 

The  above  Taylor  series  expansions  about  time  t  negiecta  alll  terms .  of 
order  two  and  higher.   The  changes  in  position,  velocity,  and  radius 


are  the  result  of  nonzero  velocity  v  e  dr/dt,.  acceleration  a  E-dv/dt, 
and  radius  change  rate  v  =  dji/dt,.  respectively.. 

Since  the  transport  equation  is  a  balance  equation  the  number,  of" 
bubbles  in  a  volume  element  d^x ''  =  d^r''d  v'dJ?- '   at  time:  t^-t- dt-,  .must 
equal  the  number  of  bubbles  that  were  in.  the  volume  d  t  =;  d  rd  vd5, 
at  time  t  plus  or  minus  the  bubbles  that  were  generated. or  removed  by 
any  internal  sources  or  sinks..  Thus  an  equation,  fcir.  the.  conservation^ 
of  bubbles  is  vrritten  as : 


i;;(r' ,v',z',t+dt)d^T'  -  i(i(r  ,v^£  ,.t)d^-r  =  S:(r  ,v,jl  ,t)d,^T  (2) 

where  S(r  ,v,5,  ,t  )d'^T   is  the  rate  at  which  bubbles  were,  introduced  into 
the  volume  element  d'^x  about  (r,.v',.2,)  by  distributed  sources.. 

Expanding  i[i'  =   ij;(r'' ,v' ,5.' 'yttdt)  in  a  Taylor,  series  about .  il'.  =- tjj(r,v,Jl,t) 
and  keeping  only  terms  thro.ugh  the  first  order  yields  :• 

ip*    =   i)  +   (3ii;/3t)dt  +  (v  +  V  ij;)dt  +  (a'  -  V  i|;)d±: -h  v(dij>/3t)dt-:      (3) 

where 

->• 

V  (•)  E  &(-)/ax,  3(-)/3y;:  3(-)/3z 

-> 

V  (•)  E  3(-)/3v  ,  3(-)/3v  ,.  3(-)/3v  .. 
V  K         y  z 

The  volimie   element  d"^!*'  was   transformed   into  terms   of  d^x  by  using 

a  Jacobian.      The  results   of  this   transformation  yielded. 

d7^r    ^    3(r%v',r')    ^7^  (4) 

3(r,v»il) 


with 


3(r',v',£    )    ^    1  +   (v    ••a)dt  t   (3v/3£)dt...  (5) 

3(r,v,£) 


Substituting  (3),  (M-),  and  (5)  into  (2),.  dividing  by  d"^Tdt 
and  taking  the  limit  dt  -»■  0  yields  the  bubble  transport  equati-on 
[Ref.  5]: 

It  should  be  noted  that  this  transport  equation  is:  di-ififerenir-  from 
the  usual  transport  equation  encountered  because  the  term  (V  -a  -h  9v/dt) 
is  generally  nonzero  as  the  bubble  radius  change  rate,.  v(r,y",£,tO  ,.  and: 
acceleration,  a(r,v,£,t),  are  functxons  of  b\±)ble:  radius  and  bubble 
velocity  respectively. 

The  distributed  sources  or  sinks  are  contained  in  the  term 
S  =  S(r,v,S,,t).   This  term  can  also  include  the  rate  at  which  biibbles 

•  •  7  7 

are   introduced   into  d  t      and  removed  from  d  t     by   scattering  of P  smallV 
scale  turbulent  eddies   or  hj  bubble-bubble   interactions..     It.  was   assumed. 
the   S  was   independent  of  i|j  or  proportional  to   ^..       This,  simplication 
yielded  a  first   order,   partial  differential  bubble  transportr  equation.. 


ni.      FORMAL   SOLUTTOIJ  TO  THE  TTTAIISTOFT  EQUATION 

The  bxibble  transpart  Eq^..    (6)    i'n  the  forra  ofr  a^  firsir  onder/ partial 
differential  equation  was   salved  using  the  method;  of;  characteristics 
[Ref.    51.      This   was   equivalent  to  repladhg  the  transporf- equation 
with   a  set  of  s imultaneous  ^  first  order,,  ordinary  differential  equations- 

The  characteristic  equations   were   determined-  fo. he  the  bubble 
dynamic  equations : 

dr/dt  =  V  r(t  )    =  r"    =^  (x.  ,y    ,z  )  (7a) 

O  O  o      o      o 

dv/dt  =   a  v(t    )    =  V     :^   (v      ,v      ,v      )  (7b) 

o  o  xo      yo      zo 

d£/dt   =  V  £(t    )    ^  £  (7c) 

O  Q 

dt^/dt  =  S   -   Z^^  i|j(t    )    =   ^,     ..  (7d) 

The  model   chosen  determined  the   acceleraflon  a(r,v,2,,t-)    andr.rate   of 

biibble  radius   change- v(r^v,.£  ,t)  ,,•  thus   establishing  the.  family  of;  char- - 

acteristic   curves.      This    family  was  parameterized  by  the  set   of:  initial 

conditions    (r    ,v   ,£    )  .,  and  di   (r    ,.v    ,£   )    was   the  initial  distribution., 
ooo  ^o      ooo 

The  general  form  of  the   distribution,  function  V7as.  obtained  by 
integrating  Eq.    (7d)    and  making  use   again  of  a  Jacob ian..     The  formal 
solution  for  the  bubble   distribution   is    [Ref..  5]:: 

i>iv,v,!i,t)d^z   =   ^(^Q^^o'^Q't)  d'7x^+  /   S(r^',y'^£",t")d7:,"dt" 

•^o 
The  integral  here  is  a  line  integral  along  a  characteristic;  curve  as 

are  all  the  integrals  in  this  thesis..  Equation  (8)  is  equivalent  to 

saying  that  the  number  of  bubbles  In  dJr   about  (r..,v,£)  at;  timer  t:  is. 

the  number  of  bubbles  that  were  originally  in  d'^x   plus  or.  minus  those 

added  or  subtracted  by  sources  or  sinks  while  proceeding  along  the 

characteristic  curve. 
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IV.   THE  MODEL  FOR  THE  SOLUTIOIf  OF  THE  BUBBCE  TRAIISPORT  EQUATION 

Before  solutions  to  Eq.  (6)  could  be  obtained.  expH'catly,.  a-. model. 
of  the  upper  ocean  had  to  be  specified..  This:  rao:deL  must  provide  explicit 
expressions  for  the  acceleration  a(ryV,.£  ,.t)  ,.  the  rate  of^  change,  of. 
bubble  radius  v{r,v,i,t),   the  bubble  distribution  S(r,v,£,,t)  and_the 
circulation  field  V(r,t).   A  Cartesian  coordinate  system  was _ chosen 
with  the  origin  located  at  the  surface  of  the  ocean,  and.  the  z.-axis 
vertical  and  positive  upward.   Th.us  alX  depths,  were  taken  to  be  negative. 

A.   THE  BUBBLE  RADIUS  CHAliGE  R/JE 

The  rate  of  change  of  bubble  radius  for  biibbles  circulating  in  the. 
ocean  is  caused  by  either  gas  diffusion,,  a.  change  in  the  hydrostatic 
pressure,  surface  tension,  or  a  combination  of  these  phenomena..  This 
interaction  of  effects  was  analyzed,  by  assuming  the  bubbles  to.  be: 
spherical  containing  an  ideal  gas ►   From  the  ideai  gas  law 

V=|..i3=B|T  (9) 

where  R  is  the  gas  constant,  T  is  the  temperature,.  P  is.  the  pressure, 
and  V  is  the  bubble  volume.   By  taking  the  time  derivative  of  Eq.  (9) 
and  assuming  isothermal  changes,  an  expression  for  v,.  the  change  in 
bubble  radius  was  obtained: 
d£    £ 


^  =  d¥-  3 


fl  dn    1  dP 
I  n  dt   P  dt 


(10) 


The  term r—  contains  the  gas  diffusion,  effects  whiXe  the  effects  of 

n  dt 

the   change   in  hydrostatic  pressuire   and  surface  tension  ar.e_-     contained 

m  the  term  —  -r—  . 
P  dt 


The  hydrostatic  pressure  P  is  the  sum  of  the-  atmospheric  pressure,  . 
P  ,  and  the  fluid  pressure,  -p-  gz..  The  surface  tension  pressure  is 
7-^1%    ,  where  p   is  the  density  of  sea  water,,  p   =1.03  x  ID  ^ '  Kg/meter-^ ,  . 
and  Y  is  the  surface  tension  of  sea  water,,  approximalrely -equal -to  0.7'+- 
newlrons/meter.   Therefore 

p  =  p   _  p  gz  +  ZyA  (11) 

Q      O 


and 


dP  2y  d)i 

dr=  -  ^oS^z  "F'd^ 


--^^^^-TY^-  (12) 


By  using  a  simple  model  for  gas  dxffusioTi  [Ref :. .  8 D  ttie:  following 
expression  vfas  obtained  for  dn/dt 

^=  -  5^£3  (p  _  p^)  (;L3) 

where  5  is  a  constant.   A  value  for  6  had.  been  determined  from  experi-- 
mental  data  to  be  aboui:  10~^  meters/second..  r-t_  was  assumed  -  for.  this 
model  that  the  dissolved  gas  concenlrration  gradienir  was  confined  to 
a  thin  shell  at  the  air-water  boundary  because  of  the  bubble 's: speed. - 

An  expression  for  v  was  obtained  by  substituting  Eqs.  (11} , . ( 12) , 
and  (13)  into  Eq.  (10)  [Ref.  5j.. 


=  ^  =  [^„il^  "r   36RT(z£  -  r)3/r3(D-z)ii  +  2r]  (14) 

where 


^   =   dt        --z' 


and 


r    =   2Y/p^g  -   1.1+7   X   10    5  meters^ 


D   E   P    /p    g   -   10  meters. 


If  2.  is  greater  than  30  microns  (I  micron  (ix)  =  lO""^'  meters), the 
surface  tension  terms  can  be  ignored,,  i.e..,.  Y  -   0.. 
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B.   THE  BUBBLE  EQQATIOrr  OF  MOTION  AND  ACCELERATION 

A  bubble  that  is  ia  motion  with  respect  to  a  stationary  reference 

system  will  experience  several  forces  acting  on  i"t. .  First  there,  is ' 

the  force  of  the  water  acting  on  the  bubble  to  pull  itr.  along. .  Iffthe 

bubble  were  completely  entrained  by  the  flui'd  so  that::  the  bubble's: 

velocity  (v)  were  equal  to  the.  circulation  field  velocity  (V),.the 

entrainment  force  would  be  p  ffdV/dt.  where  a   is  the:  volume  of.  the:-. 

o 

bubble.      Complete  entrainment.  can.  occnir'  only-  If;  p.:  is:  equal'  toip:  .  . 

Since  this   is  rarely  the   case^  the  bTobble.  will  raaihtain:  a.  relative- 
->•       ->-       ^ 
velocity  u  =  v  -  V  with  respect,  to  the  water..     This  relative  velocity 

-> 

produces  a  drag  force,  F  y  acting  on  the  bubble  to  slow  it  down.. 

As  the  bxibble  moves  through  the  water,,  the  biobble  pulls  vrater. 
along  with  it.   The  amount  of  wa±er  pulled,  along  with  it  is  equal^to 
some  fraction,  B,  of  the  bubble's  volume..  From  experimental  data  B  is 
usually  taken  as  0.5.   This  phenomena  is  called  the  mass  adhering  force 
and  is  equal  to  -(Bcrp  du/dt)..  The  bubble  will,  also  experience  .  a  net 
bouyant  force  v/hich  is  equal  to  cfg(p   -  p)  z,.  where  g  is  the  gravita- 
tional acceleration. 

By  using  Newton''s  second  law  these  forces  were,  combined  to:  yield 
an  expression  for  the  acceleration  a  CRef..  5]: 

^  "  op  (B+p/p  )  *   (6+p/p  )  dt  *   (B+P/P  )  8^ 
o      o         o  o 

where 

^-$.   (U   )V  ..    ■  (16) 

dt   &t       r 

In  order  to  obtain  an  explicit  expressi'on  for.  the  acceleration:  a,  . 
it  was  necessary  to  derive  an  expression  for  the  drag  force  F_^.   Various 
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models  for  the  drag  force  were  considered,  V7ith  the  rnodels  being 

classified  by  using  the  Reynold -'s  number  R  —  2£up  /n'  v^herc  n  is' 
J  °  J  e       o 

the  viscosity  of  water  ( n  ^  1Q~^  Newton— second/metex^^) . . 

For  bubbles  with  a  radius  less  than,  or  eqiial  to  lOD;  irricrons 
{I   <_   100  ),  the  ReynoWs  number  is  less  than  one  (R  <!)  andithe 
drag  is  viscous  according  to  Stokes  Law:- 

^D  ^  *Stokes  =  -'^"'="  ("' 

where  K  =  677  in  the  presence  of  surface  active- material,  and  M-'W"  in:  the 
absence  of  these  surface  active  materials.. 

The  drag  was  considered  to  be  both  viscous;  and.  turbulent- if:' 
loop  ;^il<_  1500  u  or  1  <  R  <  7QQ..  Thus  the  drag  force  was 

^D   "   ^Stokes   ^   ^Turb"  ^^^^ 

-»  ->- 

Here   F^^  ,         was   given  as   in  En..    (17)   with   5Tr   <   K   <   12Tr- and   F;^     .  ,.v/hich 
Stokes  ^  ^  —      —  Turb. 

arise   from  the   energy  dissipated  in  the  turbullence  of  the  wake,,  had  the 
form: 

F„     .     =   -0.5   K^    s,    uu  (19) 

Turb  f^o   1 

where  K^  is  a  drag  coefficient  and  s,  is  the  bubble  surface  area" 

f         '^  1 

covered  by  the  turbulent  wake  (s   =  t\1    ).      If  there  were  no  surface 

->  -> 

active  materials  present  F„  would  equal  F„   ,   ,.  with  K  ="  12it.  . 

D        ^     Stokes 

Since  the  oceans  contain  surface  active  materials ,.  the  drag-force 
contains  both  a  viscous  drag  and  turbulent  drag  term.   By  noting  that 


p/p   <<  1  and  that  a  =  ^irl^ /2   for  bubbles  with  £  <  1500  y , .  the.  acceler-- 
o 


3^ 
ation  can  be  written  as  [Ref..  5l; 


a  =  -  (4-^u  -  (-^)uu  +  (^)S-  +  (1  t  ^)dV/dt:  (20 
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with 

a  =  —  - —  -   8.75   X  10~^  me-ters/sec^ 

2  Bp^ 

and 

If  the  bubbles  were  generated  by  a  surface  SDurce-  such  as  raifi,  . 
wave  action,  etc.,  the  bubble  radi'i  are  generally  less"  than  150:  p 
[Refs.  3  and  9l  and  the  drag  farce  is  pred.Dra±hat"ely  viscous. .  Ihr^this' 
case  C  should  be  set  equal  to  2rera.. 

If  the  bubbles  originate  at  a  bottom  soiurce,.  the  radii,  are  ih' the 
range  of  500 y to  700 u  [Ref.  ^J   and  the  drag  should.be  considered  to  be 
predominately  turbulent  with,  a  set:  equal  to  zero.. 

C.   BUBBLE  SOURCES 

The  bubbles  in  the  ocean  might  physically  be  considered  to  have 
originated  from  rhree  t.5rpes  of  sources..  The  first  type  included. 
bubbles  that  originated  near  the  surface  and.  were  caused,  by  v/ave. 
action,  rain  or  snow.   The  distribution  of  bubbles  originating' from 
these  phenomena  was  analyzed  in  Ref .  10 .. 

The  second  type  of  bubbles  are  those  nucleated  by  distributed 
particles  of  material  as  well  as  bubbles  that  originated,  from  biological 
activity  such  as  photosythesis  or  bacteria  and  are  distributed  through- 
out the  medium. 

The  third  type  of  bubbles  includes  those,  biibbles  generated  on  the 
floor  of  the  ocean.   These  bubbles,  are  often  sediment— initiated  bubbles  : 
and  measurements  indicate  these  bxibbles  to  be  beti-jeen.  M-50  y.-- 800  y  in 
radius  with  a  peak  around  550  p  [Ref ..  4j..   It  was  this  category  of 
bubbles  for  which  the  distribution  function  was  evaluated.   It  was 

13 


considered  to  be  mathematically  convenient  to  introduce  the  biibbles 

via  boundary  conditions  rather  than  to  include  them  m  S(r  ,v,Jl  ,t) . . 

->■  -> 

t|^   (r  ,v  ,£  )  represented  the  bubble  distributxcn  originating  in" 

the  bottom  sediment. 

The  model  for  the  circulation  field  was  chosen"  to  be  s   quiescent": 
fluid  dictating  that  V  =  0.   For  this  model  the- bubbles  could.:  rise. 
only  via  buoyancy  and  the  characteristic  curves  are  verticals  upward . . 
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V.   RESULTS 


For  this  thesis  the  buLble  transport:  theory  was  app^lied.  to^  a: 
one- dimensional  quiescent,  steady  state  acean  with'  a  hottom  source 
at  197  meters.   The  solution  and  appraxiinate^  solutionrs.  are  givenr. . 
The  characteristic  equations  were  integrated  numerically  to  test' 
the  accuracy  of  the  approximate  analytical  solutions^.  The  calculated 
bxibble  distributions  were  compared  with  experimental,  data. . 

In  this  steady  state  model  the  velocities  ar^e  vertical  upward,. 
and  all  functions  depend  spatically  on.  depth  only..  The- bubble  transport" 
equation  (Eq.  6)  can  be  written,  as:: 

v|i+   a|^t   v|f  =   -r  ^     -HS..  (21) 

dz  dv  dJl  t^ 

z 

Here  ({)(z,v,£)  is  the  distribution  function,.  S  is-  the.  distributed  source  . 

function,  v  is  the  vertical  velocity,,  v  is.  the  rate,  ot  change:  of:' bubble 

radius,  and  S   =  8a/ 3v  +   8v/a£.. 

If  z  is  chosen  as  the  independent  parameter,,  the.  characteristic 

equations  are  [Ref.  51: 

^  =  -=  -^-T-^T-  (22a) 

dz   V     £    ^     Sv 

dZ        V    £2  +  35RT(zS,  -  r)/v 


dz   V     3(D-z)£  +  2r 


(22b) 


and 


^=   --     l^^      t^S..  (22c) 

dz     V   t     V 

Since  the  characteristic  curves  are  vertically  up>;ar-d:  along;  ttie:  z- axis ,  . 

the  characteristics  or  boundary  conditions  are  v  -  v  ,.  Z   =   i    ..and 

^  o       o 

<|»  =  (i>(z    ,v  ,£  )  at  z  =  z   <  0. 
o  o  o  o  — 
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Equation  (22c)  was  iategrated  to  yield 


z 
(})(z,v,5.)  =  (j3(z  .^qT^q)  expf-  /   i:^(z',.v'-,Jl')dz'/v'_[ 

/z   ^o 
S^(z",v'\;i")dz"/v"  J  S(z' ,v' ,Ji')d2J/v    (23) 


with  the  integration  taking  place  along  the  char  act  eris'tic  curves: 
{z,v(z:z  ,v  ,£  ),£(z:z  ,v  -£  );:  z  <  0 } ..  By  using"  the  Jacohian .  Eq. .  ( 23) 

OOO         OOQ       — 

was  rewritten  as  [Ref .  5] r 

V(i(z,v,£)dvd2,    =   V  (b(z    y\r   ^ji    )dv-  dfi, 

^'  Q^OOOOQ 

+    /     S(z',v',JiMdz'dv'dil''..  (2^) 

yzo 

Equation  (24)  is  interpreted  as  a  canservati'oTi  equation  which  states 
that  the  flux  of  bubbles,  V(^ ,   at  the  point  (z,.v,Jl)  is  equal,  to  the 
flux  of  bubbles,  v  A  ,  at  the  point:  (z  ,.v  ,.£  )  plus  (or  minus-)  any 

O^O  ^  0   0   0*^ 

bvibbles  added  (or  substracted)  by  distributed  aour.ces  along  the 
characteristic  curve. 

Prior  to  the  calculation  of  the  distribution  function  in  either 
Eq.  (23)  or  Eq.  (24)  it  was  necessary  to  establish  a  mapping  between 
points  (z,v,£)  and  (z'jv'ji*")  on  the  characteristic  curve-.  This  was 
done  by  integrating  simultaneously  the  characteristic  equations  (Eq. 
22)  to  obtain  v(z,£;  z%v',£')  and  £(  z,£;z' ,v' ,£ '  )  for  various. 
z' ,v' ,  and  £'  and  substituting  these  relationships  into  Eq.  (23).. 
Since  a/£^  >>  Ij  it  was  very  time  consuming  to  integrate  Eq..  (22a-) 
numerically.   Referring  to  Fig.  1,.  it  is  seen,  that  an:  appropri-ate 
approximation  is  to  equate  the  velocity  and  the:  terraihal'  velocity. . 
This  permits  semi-analytical  solutions.. 
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This   approximation  can  be   verified  by  rewriting  Eq.    (22a)    as 
[Ref.    5] 

^  ^  _    (1.)    (v  _  V    )    (v   -h  v'     ) 


where 


and 


^T   = 


-1  +   Cl  +  ^(v^]_/v-T.2^^ 

2(v      /v      2.) 
^    Tl^    T2    ^ 


1   t  Cl  +   ^^^Tl/^T2^^^ 


I  y^ 


y/2 


(25) 


(25) 


(27) 


'Tl 


3   a 
t  i. 


(28a) 


is   the  terminal  velocity  in  the   case  where  viscous   drag  dominates 
(£    <_  100    ti). 


T2 


£1 
0  ? 


(28b) 


is  the  terminal  velocity  in  the  case  where  turbulent  drag:  dominates 
il   >   M-00  y).   Thus 


^T  = 


V    for    £  <_  IDO 


V    for    I   ^  il00 


(29) 


Equation  (25)  was  integrated  using  fourth  rank  optional  Runge  Kutta 

[Ref.  11].   These  numerical  results  are  presented  graphically  in  Fig.  1. . 

It  is  observed  that  when  a  bubble  is  released-  from  a  depth  of  20.00 

meters  with  zero  velocity,  it  ^-jill  reach  terminal  velocity  after,  having 

traveled  a  distance  of  only  several^  bubble  radii'..  Thus  the  approximation 

that  V  =  V  is  considered  to  be  a  valid  approximation  for  this  model. 

This  approximation  is  convenient  as  Eq.  (22a)  is  very  difficult  to 

integrate  numerically. 
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By  using  the  approximation  v  =  v;^^  the  bubble  dislmibution  function 
can  be  written  approximately  as  r. 


where 


<^(z,v,il>  =  5(v  -  v-^)$(z,l) 


'(z,Jl)  =   /  ({)(z^v,£)d^r 


(30) 


(31) 


is  the  "depth-radius"  distributLon  function..  By  s-ubs-ti±utlhg- Elq. .  (31) 
into  Eq.  (24)  and  integrating  over  v,.  an  expression  for^  $  was:  determined 
as  [Ref.  5] 


To 


^(Z,£)  =  $(z  ,1     ) 


where 


and 


S(z',il')  E 


dz' 


(32) 


S(z%v%ii')dv'' 


A^]'- 


dl      3 


etc.. 


(33a) 


(33b) 


In  order  to  calculate  $(zy5.)^  it  was  necessary  to  determine  8£'/'3£ 
and  il'(z';z,Jl)  either  analytically  or  numerically.   Notre  thatr 
9£'/35,  =  l/(3£/32,'),  and  that  the  approximate  differential  equation 
that  9V9£  '  satisfies  is  obtained  by  setting  v  =  v   in  Eq..  (22b)  and 
taking  the  partial  derivative  with,  respect  to  Jl'CRef..  5.]: 

r 


9_     92_ 

37    aT* 


mvi'^t-)]/[ 


3(D-z) 


2r 


I  J 


[-( 


35  RT 


(t 


•^      dX 


3(D-z)    -h^- 


(34a) 


where 


-1/2 


3v^_    3v^ 
3£  £ 


1   +    4 


V       v2 

Tl 


T2 


T 

"IT 


(34b) 


Furthermore,  with  the  approximation  v  =  v  ,  then  ^=   ii„,  and  Eq.  (22b) 
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becomes 

dJL     l^   +  36RT(z£  -  r)/v 

dz         3(D-z)£  +  2r  ''''^'' 

Equations  (3M-)  and  (35)  can.  be  integrated  mimerically  using  optimal. 

fourth  rank  Range  Kutta  CRef.  Il3 ..  These  results  can  be  substituted 

into  Eq.  (32)  to  obtain  $  =  <5   for  various  soxircBS^..  Sinrce  v  =-v_  is 

a  good  approximation,  it  is  expected  that  the  approximation.5-==£   will. 

also  be  very  good. 

To  obtain  approximate  analyiiical  results  for  comparison  -  to  the 

more  exact  numerical  results  outlined  above,,  the  approximatlbns"  of:' 

Eq.  (29)  can  be  used.   Neglecting  surface  tension,.  T  —  o  is  also  a 

valid  approximation  for  I   >    30  u..  Then  Eq..  (35")  can  be  integrated 


analytically  [Ref.  5l.   For  v  -  v   and  T  -  0 


;z.,..>=[| 


^1/3 


Tl   '   '      |D-z        2  (g/3)    D-z_ 


(36a) 


For  V  =  V   and  F  =  0 


.=.„<..-..- ((^)"'-""--SfeT7ir»-.te 


z'   1 


((g/B)0^/^   ^ 
_3Dfl-f?^V/^ll]  ''\.  (37a) 


(.-(] 


The  inverse  transformation  are 


i  =  ^^^(z';z,<i)      •  (36b) 


and 


I   =   S.^^iz'',z,l)  ■  (37b) 

re  spe  ct 1 ve ly .   Then 

3£'    ,D-z  .   £^  f^f-^^ 
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and 

3£»    _  /D-z    I        /   I     \ 

"31"  '  V^^^l        Ur-]  (37c) 


respectively.        If  gas   diffusion   can  be  negiectred,.  6    —  0,then    E'qs..(35) 
and   (37)   reduce  to: 

1/3 
^=     \l^-r-\  (38a) 


1/3 
'''       '^-^     ^  .  •  (88b) 


3i2. 

Figures  (2)-,  (3),  (*+)»  (5),  and  (5)  demonstrate  the-:  results. obtained 
by  integrating  Eq.  (35)  numericalXy  and  coinparing^  £  to' tte  analytical- 
approximations  of  Eqs-  (36),  (37)^  and_(38)..  Figures  (2.)  , .  (  3) ,  .  and  (4) 
were  draiNm  for  a  depth  of  20  meters  and  Figs..  (5)  and  (5)  were  drawn 
for  a  depth  of  197  meters.   Figures  (2)  and.  (5^)  demons1:rate- that .  £  . - 
is  a  good  approximation  for  I     when  £.  <  100  u...  However, .  this  approxi-- 
mation  breaks  down  when  Jl  <  30  \i:.     This  is  because  the  effects  off 
surface  tension,  which  are  important  in  bubbles  of.  £  <  30  y.,  .are; 
neglected.   Figures  (3),  (4),  and  (6)  demonstrate  that  £   isagood.. 
approximation  for  2.  when  I   >_   UOO  \i..      From  Fig.,  (5),.  it  is.  observed 
that  I        is  not  a  good  approximation  for  £..  This  is  because  £ 
neglects  the  effects  of  gas  diffusion.. 

The  survival  threshold  for  bubbles  is  differenf  depending  on 

whether  the  bubbles  were  released  at  a  depth  20  meters  or  197  meters. 

In  the  20  meter  case,  bubbles  of  £  >  100  \i   will  reach  the.  surface,. 

while  in  the  197  meter  case,  bubbles  of  £>  LOOO  meters,  will,  surface.  . 

This  difference  is  caused  hy   gas  diffusion  being  dominate,  over,  the . 

decrease  in  hydrostatic  pressure   Since  gas  diffusion  is  proportional' 

to  (P-P  )/P,  it  will  be  constant  for  z  <  30  meters  and  will  only 
o  —  -^ 
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decrease  for  z  _>  30  meters.   Thus  bubbles  released  at  19T  meters  are 

exposed  to  a  high  gas  diffusion  rate  for  practically  all  their-  travel,. 

while  bubbles  released  at  20  meters  are  exposed  to:  a  decreasing  gas 

diffusion  rate. 

Finally  the  bubble  distribution  function  was  calculated  for.  each 

of  these  approximate  solutions  to  Eqs..  (22a)  and  (22b)j  (3M-)  and;  (3  5),. 

(36),  (37),  and  (38).   It  was  convenient  to  introduce  the  initial- 

distribution  function  as  the  boundary  condxtlon.  0(z;z:  ,£  )  instead- 

"^  o   o: 

of  as  a  distributed  source.   Since  s(z'',X'')  is  zero,,  the  solution  can 
be  written  as 

V 

$(z;z  ,£  )  =  ^  <l(r  ,£  )  l/0£/3£')..  (39a) 

o  o    V       o   o 

For  V  =   V  ,  Eqs.  (34)  and  (35)  must  be  integrated  numerically  to. yield 

<I>^(z'z  ,il  ).   For  V  =v^,(£  <  lOOy)  and  r  -  0 
T.   '  o   o  Tl   — 

$  =  <&,Cz;z  ,£  )  =   l^f^]   'P(z  ,1  );  (39b) 

1  o   o     \D-z  J  o   o 

for  V  =^2^^  -  ^^°^^ 

$  =  $^(z;z  ,1    )  =f^-^-     $  (z  ,X  );:  (39c) 

2  '  o  o    I  D-z  i        o  o 

and  if  the  effects  of  gas  diffusion  are  not  considered 

V  ;     \  1/3 

$  =  $  (z;Z  ,i    )    =      -^$(z  ,£  )  {^^]  (39d) 

3  o   o     V      o   o   VD-z  / 

Here  V   is  computed  using 

in  Eq.  (26). 

The  above  distributions  were  calculated  using  an  initial  distribution 

function,  ^(z  ,£  ),  obtained  from  Ref..   4  ..  The  authors  of  Ref;  4 
o  o 
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determined  the  distribution  of  bubbles  rising  from  "the  highly  organic,. 
anaerobic  sediment  on  the  bottom"  of  Saanich  Inlet:,.  British  Columbia. - 
The  depth  of  the  water  was  197  meters..  The  density  of  these  rasing 
bubbles  was  determined  by  using  stationary  echo-soundeTS  operating 
at  various  frequencies  in  the  kilohertz  range..  Figure  7a  is  a-  histogram- 
that  presents  the  presumed  bubble  distribution,,  ^(z    ,1    )  at  a.  depth  of~ 
197  meters  as  presented  in  Ref.  M-,   Figures  (7b)  through  (7g)  demon — 
strate  how  this  initial  distribution  changesr  as  the  bubbles  rise  in." 
height.   It  shoiild  be  noted  that  height  is  measured  from  the  ocean 
bottom  upward  to  the  surface. 

The  different  histograms  at  each  height  above  the  ocean  floor  shows 
how  the  bubble  distribution  is  changing  depending  on  the.  models  chosen 
to  be  used  for  the  calculation  of  £  ..   In  the  model  vrhere  gas  diffusion 
as  well  as  hydrostatic  pressure  changes  are  considered^,  the  bubble 

distribution  <5(z:z  ,£  )  shifts  toward  smaller  size  bubbles  as  heightr, 

o  o  ° 

increases.   It  can  also  be  observed  that  the  width  ofT  the  intervals 

is  increasing  with  height.   Finally  all  the  bubbles  disappear  entirely 

at  a  height  of  about  100  meters. 

In  the  small  bubble  approximation  $  =-  <J>  ( z  ;:z  ,X  ),.  where  the  drag 

force  is  considered  to  be  viscous  and  v  =  '^m-i  »■  "the  bubble  distribution 

appears  to  spread  out  by  shifting  toward  bubbles  of  both  smaller  and 

larger  radii.   Bubbles  with  initial  radius  less  than  900  y  continually 

decrease  in  radius,  while  bubbles  with  an  initial  radius  greater  than 

900  \i   vrill  initially  decrease  in  size  but  as  the  effects  of  gas 

diffusion  reduce,  these  bubbles  expand..   In  the  turbulent  drag^model,. 

v  =  v^^ ,  $  =  $-(z:z  ,£  )  shifts  toward  smaller  radius  bubbles  because 
i  z  I  o   o 

of  the  effect  of  gas  diffusion  and  disappears  at  a  height  of  about  125 
meters . 
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When  gas  diffusion  is  neglected,  the  bubble  distribution 

4>  a  ^^(ziz  ,1    )    shifts  tovj-ard  babbles  of  larger  radii'..  This  v/as  to 
3    o  o  ° 

be  expected  because  as  hydrostatic  pressure  decreas-ed,.  the  bubbles 
would  increase  in  size.   The  bubble  distribution  calculated,  for. 
$  (z;z  ,Jl  )  compares  quite  well  with  the  experimental  data,  ofr 
McCartney  and  Bary  in  Ref .  M-.   It  appears  as  though  gas  dn  f  fusion 
was  somehow  not  important  in  the  bubbles  they  measiored..  These" 
bubbles  originated  at  the  bottom,  so  perhaps  they  had  some,  k-ind.  of. 
coating  that  impeded  gas  diffusion  across  their  surfaces.. 
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VI.   CONCLUSION 

Calculating  bubble  distributions  in  the  ocean  by^  aolving-  the: 
bubble  transport  equation  appeared  to  be  quite  feasible..  Altho^lgh 
the  solution  of  the  transport  equation  can  be  q^lite  involved^,  a  careful' 
choice  of  approximations  often  permits  approximate  analytical  solutions. 
Otherwise,  the  equations  can  usually  be  numerically  integrated^.  Eor 
the  model  chosen,  analytical  expressions  were  obtained  for  several 
approximations.   All  of  these  expressions  yielded  results  that:  were 
in  general  agreement  with  experimental  data  or  with  iniruitive  reasoning.. 
It  was  seen  that  the  drag  force  model  chosen  and  the  inclusion  of  g"as 
diffusion  greatly  influenced  the  results  obtained.. 

The  assumption  that  the  velocity  of  the  bubble  was  equal  to.  the 
terminal  velocity  should  be  quite  valid  as  the  data  presenired.  on  Fig.. 
1  demonstrated.   V/ith  this  approximation,  one  can  easily  calculate  the 
bubble  radius  at  any  depth.   The  radius  is  strongly  dependent  on  the 
model  chosen  for  the  drag  force,  gas  diffusion  and  surface  tension.. 
The  choice  of  the  drag  model  determines  the  equation  used  to  calculate 
the  terminal  velocity. 

Finally  the  distribution  function  for  a  bottom  source  was  calculated. 
The  data  generated  by  this  theoretical  approach  agrees  well  vjith  the 
limited  experimental  data  available.   The  distribution  calculated  for 
the  no  gas  diffusion  model  yields  data  that  is  in  gcod  agreement  v;ith 
the  experimental  data  in  Ref.  ^■. 

From  the  limited  amount  of  experimental  data  available,,  it  appears" 
that  gas  diffusion  is  unimportant  in  bubbles  originating  in  bottom 
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sediment  and  transported  to  the  surface  >   It  is  felt  that:  the  .type  of 
gas  in  the  bubbles  would  greatly  influence  the  amount  of  diffusion  that' 
would  take  place  [Ref.  12].   Since  a  fluid,  super-sa1:ura1:ed  v/ith  the 
same  gas  that  is  within  the  bubble,  is  the  only  fluid  in  which  gas 
diffusion  is  not  considered  the  dominating  effect,,  there  must  be.  some 
mechanism  inhibiting  gas  diffusion  from  sediment  bubbles..  Further 
investigation  into  this  phenomena  seems  warranted.. 
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1.0  2.0  3.0  4.0 

HEIGHT  OFF  BOTTOM  IN  METERS 

CHANGE  IN  BUBBLE  RADIUS  AS  A  FUNCTION  OF  DISTANCE 
BUBBLE    HAS  RISEN  FROM    Zo=-l97  METERS  A  COMPARISON 
OF  APPROXIMATE  ANALYTICAL  SOLUTION  TO  EXACT 
NUMERICAL    SOLUTION. 

Figure    5 
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CHANGE  IN  BUBBLE  RADIUS  AS  A  FUNCTION  OF  DISTANCE 
BUBBLE  HAS  RISEN    FROM    Zo  =-197  METERS  A  COMPARISON 
OF  APPROXIMATE   ANALYTICAL  SOLUTION  TO    EXACT 
NUMERICAL  SOLUTION. 

Ficrure   6 
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DISTRIBUTION  OF  BUBBLES   ORIGINATING  IN 
SEDIMENT  AT  DEPTH    197  METERS. 
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Figure   7a 
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COMPUTER  OUTPUT 


HEIGHT  OFF  THE  BOTTOM  = 


ZS.OOQ    METFTIS 


XLO 

RADIUS 

PHI 

RLT3 

FH^I3: 

400.000 

0.0 

0.0 

417.53  5 

O'.O 

465.000 

240.062 

0.924 

485.385 

0..9  29 

475.000 

253.600 

0.925 

495.823 

a..929 

485.000 

276.132 

0.926 

5  0  6.261 

a..930 

490.000 

284.593 

6.482 

511.480 

6..509 

500.000 

301.010 

6.436 

521.919 

6.510 

510. 000 

316.856 

6.490 

532.357 

6.512 

515.000 

324.596 

31.530 

537.576 

31.-633 

530.000 

347.201 

31.552 

553^234 

31.643 

545.000 

369.032 

31.572 

56  8.892 

31.653 

555.000 

383.230 

15.792 

579.330 

15r.83  0 

575.000 

410.924 

15.802 

600.207 

15..835 

595.000 

437.342 

15.812 

621.083 

15.840 

600.000 

444.469 

3.721 

626.303 

3.727 

635.000 

489.908 

3.724 

662.837 

3.729 

660.000 

521. 523 

3.726 

688.933 

3. .730 

670.000 

534.007 

0.932 

699.371 

0.933 

705.000 

577. 089 

0.932 

735. 9J6 

0,933 

730.000 

607.364 

0.933 

762.002 

0.933 

740.000 

619.375 

0.933 

772.440 

0.933 

775.000 

661.022 

0.933 

803.974 

0.,.934 

800.000 

690.444 

Q.933 

835.070 

0.934 

810.000 

702.146 

Q.G 

845.509 

Qu-0 

XLO 

RLTl 

PHIl 

RLT2 

PHLr 

400.000 

338.683 

0.0 

206.200 

0-.,0 

465.000 

431.628 

0.879 

297.174 

0..938 

475.000 

444.714 

0.879 

310.439 

0.938 

485.000 

457.588 

0.879 

323.562 

0.-938 

490.000 

463.954 

6.155 

330.075 

6.564 

500.000 

476.556 

6.155 

34  3.009 

6.564 

510.000 

439.000 

6.155 

355.830 

6.564 

515.000 

495. 169 

29.894 

362.201 

31.881 

530.000 

513.480 

29.894 

331.168 

31.881 

545.000 

531. 533 

29.894 

399.935 

31.881 

555.000 

543.442 

14.947 

412.346 

15.940 

575.000 

566.997 

14.947 

436.950 

15.940 

595.000 

590.247 

14.947 

461.293 

15.940 

600. COO 

596.018 

3.517 

467.342 

3.751 

635.000 

636.015 

3.517 

509.318 

3.751 

660.000 

664.219 

3.517 

538.952 

3.751 

670.000 

67  5.43  0 

0.379 

550.735 

0.938 

7  05.0  00 

714.391 

0.879 

591.686 

0.938 

730.000 

741.996 

0.879 

620.693 

0.938 

740.000 

752.992 

0.879 

632.244 

a.938 

775.000 

791.303 

0.879 

672.467 

0.938 

800.000 

818.520 

0.379 

701.018 

0.938 

810.000 

829.376 

Q.O 

712.401 

Q..0 

38 


HEIGHT    OFF    THE    BOTTOM    = 


50.000    METERS 


XLO 

RADIUS 

PHI 

RLT3 

PHI3 

400.000 

0.0 

0.0 

438.615 

0.0 

465.000 

0.0 

0.0 

509. 890 

0.855 

475.000 

0.0 

0.0 

520.855 

0.855 

485.000 

0.0 

0.0 

531. 821 

0.856 

490.000 

0.0 

0.0 

537.303 

5.991 

500.000 

0.0 

0.0 

548.269 

5.994 

510.000 

0.  0 

0.0 

559.234 

5.997 

515.000 

0.0 

0.0 

564.717 

29. 136 

530.000 

0.0 

0.0 

581.165 

29.155 

545.000 

0.0 

0.0 

597.613 

29.174 

555.000 

0.0 

0.0 

608.578 

14.593 

575.000 

92.385 

14.316 

630.509 

14.603 

595.000 

188.945 

14.403 

652.440 

14.613 

600.000 

204.448 

3.393 

657.923 

3.439 

635.000 

290.264 

3.410 

696.301 

3.443 

660.000 

340.102 

3.419 

723.715 

3.445 

670.000 

358.626 

0.855 

734.680 

0.861 

705.000 

419.267 

0.857 

773.059 

0.862 

730.000 

459.715 

0.858 

800.473 

0.863 

740.000 

475.396 

0.859 

811.438 

0.863 

775.000 

528. 524 

0.860 

849.817 

0.863 

800.000 

565.120 

0.361 

877.230 

0.864 

810.000 

579.499 

0.0 

888.196 

0.0 

XLO 

RLTl 

PHU 

RLT2 

PHI2 

400.000 

222.533 

0.0 

0.0 

0.0 

465.000 

389.746 

0.758 

0.0 

0.0 

475.000 

408.047 

0.758 

48.238 

0.871 

485.000 

425.532 

0.758 

79.764 

0.871 

490.000 

434.014 

5.309 

93.391 

6.096 

500.000 

450.529 

5.309 

118.240 

6.096 

510. 00  0 

466.521 

5.3  09 

140.916 

6.096 

515.000 

474.345 

25.787 

151.662 

29.610 

530.000 

497.222 

25.787 

182.147 

29.610 

545.000 

519. 336 

25.787 

210.647 

29.610 

555.000 

533.722 

12.894 

228.815 

14.805 

575.000 

561.775 

12.894 

263.610 

14.805 

595.000 

589. 029 

12.894 

296.797 

14.805 

600.000 

595.736 

3.034 

304. 888 

3.484 

635.000 

641.693 

3.034 

359.687 

3.484 

660.000 

673.641 

3.034 

397.238 

3.484 

670.000 

686.252 

0.758 

411.959 

0.  871 

705.000 

729.760 

0.758 

462.363 

0.871 

730.000 

760.329 

0.758 

497.464 

0.871 

740.000 

772.455 

0.758 

511.326 

0.871 

775.000 

814. 509 

0.758 

559.139 

0.871 

800.000 

844.226 

0.758 

592.707 

0.871 

810.000 

856.048 

0.0 

606.015 

0.0 
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HEIGHT 

OFF  THE 

BOTTOM  = 

75.000 

METERS 

XLO 

RADIUS 

PHI 

RLT3 

PHIS 

400.000 

0.0 

0.0 

464.721 

0.0 

465.000 

0.0 

0.0 

540.238 

0.775 

475.000 

0.0 

0.0 

551.856 

0.776 

485.000 

0.0 

0.0 

563.474 

0.777 

490.000 

0.0 

0.0 

569.283 

5.440 

500.000 

0.  0 

0.0 

580.901 

5.444 

510.000 

0.0 

0.0 

592. 519 

5.449 

515.000 

0.0 

0.0 

598.328 

26.474 

530.000 

0.0 

0.0 

615.755 

26.502 

545.000 

0.0 

0.0 

633.182 

26. 528 

555.000 

0.0 

0.0 

644.800 

13.272 

575.000 

0.0 

0.0 

668.036 

13.288 

595.000 

0.0 

0.0 

691.272 

13.302 

600.000 

0.0 

0.0 

697.081 

3.131 

635.000 

0.0 

0.0 

737.744 

3.136 

660.000 

0.0 

0.0 

766.789 

3.139 

670. 000 

0.0 

0.0 

778.407 

0.785 

705.000 

148.353 

0.766 

819.070 

0.786 

730.000 

244.821 

0.772 

848.115 

0.787 

740. GOO 

273.440 

0.774 

859.734 

0.787 

775.000 

358.003 

0.778 

900.397 

0.788 

800.000 

410.244 

0.780 

929.442 

0.788 

810.000 

429.964 

0.0 

941.060 

0.0 

XLO 

RLTl 

PHIl 

RLT2 

PHI2 

400. COO 

a.  0 

0.0 

0.0 

0.0 

465.000 

331. 106 

0.638 

0.0 

0.0 

475.000 

360.093 

0.638 

0.0 

0.0 

485. OCO 

386, 043 

0.638 

0.0 

0.0 

490. COO 

398.163 

4.464 

0.0 

0.0 

500.000 

421.059 

4.464 

0.0 

0.0 

510.000 

442.507 

4.464 

0.0 

0.0 

515.000 

452.784 

21.681 

0.0 

0.0 

530. 000 

482.162 

21.681 

0.0 

0.0 

545.000 

509.784 

21.681 

0.0 

0.0 

555.000 

527.421 

10.841 

0.0 

0.0 

575.000 

561.206 

10.841 

0.0 

0.0 

595.000 

593.406 

10.841 

12.645 

13.575 

600.000 

601.252 

2.551 

42.242 

3.194 

635. GOO 

654.335 

2.551 

153.845 

3.194 

660.000 

69  0.674 

2.551 

212.749 

3.  194 

670.000 

704.917 

0.638 

234.349 

0.799 

705.000 

753.684 

0.638 

304.416 

0.799 

730.000 

787.663 

0.638 

350.758 

0.799 

740.000 

801.088 

0.638 

368.659 

0.799 

775.000 

847.439 

0.638 

429.073 

0.799 

800.000 

880.028 

0.638 

470.501 

0.799 

810.000 

892.959 

0.0 

486.742 

0.0 
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HEIGHT    OFF    THE    BOTTOM    = 


100.000    METERS 


XLO 

RADIUS 

PHI 

RLT3 

PHI  3 

400. COO 

0.0 

0.0 

498.412 

0.0 

465.000 

0.0 

0.0 

579.404 

0.690 

475.000 

0.0 

CO 

591. 865 

0.691 

485.000 

0.0 

0.0 

604.325 

0.692 

490.000 

0.0 

0.0 

610.555 

4.846 

5  00.0  00 

0.0 

0.0 

623.016 

4.851 

510.000 

0.0 

0.0 

635.476 

4.856 

515.000 

0.  0 

0.0 

641.706 

23.599 

530.000 

0.0 

0.0 

660.396 

23.634 

545.000 

0.0 

0.0 

679.087 

23.666 

555.000 

0.0 

0.0 

691.547 

11.843 

575.000 

0.0 

0.0 

716.468 

11.862 

595.000 

0.0 

0.0 

741.388 

11.880 

600.000 

0.0 

0.0 

747.619 

2.796 

635.000 

0.0 

0.0 

791.230 

2.803 

660.000 

0.0 

0.0 

822.380 

2.807 

670.000 

0.0 

CO 

834.841 

0.702 

705.000 

0.0 

0.0 

878.452 

0.703 

730. 000 

0.0 

0.0 

909.603 

0.7  04 

740.000 

0.0 

0.0 

922.063 

0.704 

775.000 

0.0 

0.0 

965.674 

0.705 

800, 000 

136.578 

0.680 

996.825 

0.706 

810.000 

194.535 

0.0 

1009.285 

0.0 

XLO 

RLTl 

PHII 

RLT2 

PHI2 

400.000 

0.0 

0.0 

0.0 

0.0 

465.000 

219.089 

0.517 

0.0 

0.0 

475.000 

285.770 

0.517 

0.0 

0.0 

485.000 

332.348 

0.517 

0.0 

0.0 

490.000 

351.978 

3.618 

0.0 

0.0 

500.000 

386.703 

3.618 

0.0 

CO 

510.000 

417.224 

3.618 

0.0 

0.0 

515.000 

431.336 

17.575 

0.0 

0.0 

530.000 

470.296 

17.575 

0.0 

CO 

545.000 

505.502 

17.575 

0.0 

0.0 

555.000 

527.438 

8.787 

0.0 

0.0 

575.000 

568.553 

8.787 

0.0 

CO 

595.000 

606.878 

8.787 

0.0 

0.0 

600.000 

616.113 

2.068 

0.0 

0.0 

635.000 

677.771 

2.068 

0.0 

0.0 

660.000 

719.318 

2.068 

0.0 

0.0 

670.000 

735.487 

0.517 

0,0 

0.0 

705.000 

790.452 

0.517 

54.319 

0.719 

730.000 

828.450 

0.517 

147.833 

0.719 

740.000 

843.406 

0.517 

177.412 

0.719 

775.000 

894.835 

0.517 

267.406 

0.719 

800.000 

930.829 

0.517 

324.165 

0.719 

810.000 

945. 079 

0.0 

345.729 

0.0 
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HEIGHT  OFF  THE  BOTTOM  = 


125,000  METERS 


XLO 

RADIUS 

PHI 

RLT3 

PHI3 

400.000 

0.0 

0.0 

544.643 

0.0 

465. 000 

0.0 

0.0 

633.148 

0.596 

475.000 

0.0 

0.0 

646-764 

0.597 

485.000 

0.0 

0.0 

660.380 

0.598 

490.000 

0.  0 

0.0 

667.188 

4.191 

500.000 

0.0 

0.0 

680.804 

4.  197 

510.000 

0.0 

0.0 

694.420 

4.203 

515.000 

0.0 

0.0 

701.228 

20.428 

530.000 

0.0 

0.0 

721.652 

20.467 

545.000 

0.0 

0.0 

742.076 

20.504 

555.000 

0.0 

0.0 

755.693 

10.264 

575.000 

0.0 

0.0 

782.925 

10.286 

595.000 

0.0 

0.0 

810.157 

10.306 

600.000 

0.0 

0.0 

816.965 

2.426 

635.000 

0.0 

0.0 

864.621 

2.433 

660.000 

0.0 

0.0 

898.661 

2.438 

670.000 

0.0 

0.0 

912.278 

0.610 

705.000 

0.0 

0.0 

959.934 

0.611 

730.000 

0.0 

0.0 

993.974 

0.612 

740.000 

0.0 

0.0 

1007. 590 

0.613 

775.000 

0.0 

0.0 

1055.246 

0.614 

800.000 

0.0 

0.0 

1089.287 

0.614 

810.000 

0.0 

0.0 

1102.903 

0.0 

XLO 

RLTl 

PHU 

RLT2 

PHI2 

400.000 

0,0 

0.0 

0.0 

0.0 

465.000 

0,0 

0.0 

0.0 

0.0 

475.000 

0.0 

0.0 

0.0 

0.0 

485.000 

237.512 

0.396 

0.0 

0.0 

490.000 

281.882 

2.773 

0.0 

0.0 

500.000 

344. 694 

2.773 

0,0 

0.0 

510.000 

392.070 

2.773 

0.0 

0.0 

515.000 

412.549 

13.469 

0.0 

0.0 

530.000 

466.053 

13.469 

0.0 

0.0 

545.000 

511.787 

13.469 

0.0 

0.0 

555.000 

539.423 

6.734 

0.0 

0.0 

575.000 

589.940 

6.734 

0.0 

0.0 

595.000 

635.914 

6.734 

0.0  • 

0.0 

600.000 

646. 868 

1.585 

0.0 

0,0 

635.000 

719.049 

1.585 

0.0 

0,0 

660.000 

766.966 

1.585 

0.0 

0,0 

670.000 

785.493 

0.396 

0.0 

0.0 

705. COO 

848.073 

0.396 

0,0 

0.0 

730.000 

891.037 

0.396 

0.0 

0.0 

740.000 

907. 893 

0.396 

0.0 

0.0 

775.000 

965.657 

0.396 

0.0 

0.0 

800.000 

1CK)5.929 

0.396 

114.038 

0.629 

810.000 

1021.843 

0.0 

153. 122 

0.0 
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HEIGHT  OFF  THE  BOTTOM  = 


150.000  METERS 


XLO 

RADIUS 

PHI 

RLT3 

PHI3 

400.000 

0.0 

0.0 

614.835 

0.0 

465.000 

0.0 

0.0 

714.746 

0.490 

475.000 

0.0 

0.0 

730.116 

0.491 

485.000 

0.0 

0.0 

745.487 

0.492 

490.000 

0.0 

0.0 

753.173 

3.445 

500.000 

0.0 

0.0 

768.544 

3.451 

510.000 

0.0 

0.0 

783.914 

3.458 

515.000 

0.0 

0.0 

791.600 

16.809 

530.000 

0.0 

0.0 

814.656 

16.851 

545.000 

0.0 

0.0 

837.713 

16.890 

555.000 

0.0 

0.0 

853. 083 

8.457 

575.000 

0.0 

0.0 

883.825 

8.481 

595.000 

0.0 

0.0 

914.567 

8.502 

600.000 

0.0 

0.0 

922.252 

2.002 

635.000 

0.0 

0.0 

976.050 

2.009 

660.000 

0.0 

0.0 

1014.478 

2.014 

670.000 

0.0 

0.0 

1029.848 

0.504 

705.000 

0.0 

0.0 

1083.646 

0.506 

730.000 

0.0 

0.0 

1122.074 

0.506 

740.000 

0.0 

0.0 

1137.445 

0.507 

775.000 

0.0  ■ 

0.0 

1191.243 

0.503 

800.000 

0.0 

0.0 

1229.670 

0.509 

810.000 

0.0 

0.0 

1245.041 

0.0 

XLO 

RLTl 

PHIl 

RLT2 

PHI2 

400.000 

0.0 

0.0 

0.0 

0.0 

465.000 

0.0 

0.3 

0.0 

0.0 

475.000 

0.0 

0.0 

0.0 

0.0 

485.000 

0.0 

0.0 

0.0 

0.0 

490.000 

0.0 

0.0 

3.0 

0.0 

500.000 

288.314 

1.928 

0.0 

0.0 

510.000 

372.655 

1.928 

0.0 

0.0 

515.000 

404.245 

9.362 

0.0 

0.0 

530.000 

480.123 

9.362 

0.0 

0.0 

545.000 

540.491 

9.362 

0.0 

0.0 

555.000 

575.747 

4.681 

0.0 

0.0 

575.000 

638.591 

4.681 

0.0 

0.0 

595.000 

694.510 

4.681 

0.0 

0.0 

600.000 

707.697 

1.101 

0.0 

0.0 

635.000 

793,635 

1.101 

0.0 

0.0 

660.000 

349.931 

1.101 

0.0 

0.0 

670.000 

871.653 

0.275 

0.0 

0.0 

705.000 

944.486 

0.275 

0.0 

0.0 

730. 000 

994.222 

0.275 

0.3 

0.0 

740.000 

1013.686 

0.275 

0.0 

0.0 

775.000 

1080.214 

0.275 

0.0 

0.0 

800.000 

1126.462 

0.275 

0.0 

0.0 

810.000 

1144.711 

0.0 

0.0 

0.0 

43 


HEIGHT    OFF    THE    BOTTIM    = 


ETS-OOa    METER S: 


XLO 

400.000 
465. COO 
475.000 
485.030 
490.000 
500.000 
510.000 
515.000 
530.000 
545.000 
555.000 
575.000 
595. OCO 
600.000 
635.000 
660.000 
670-000 
705.000 
730.000 
740.000 
775.000 
800.000 
810.000 


RADIUS 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0 
0.0  • 
0.0 
0.0 


PHI 
0.0 
0.0 
Q.O 
0.0 
Q'.O 
Q.O 
0.0 
0.0 
0.0 
Q.O 
0.0 
Q.0 
0.0 
Q.O 
Q.O 
0.0 
0.0 
Q.O 
Q.O 
0.0 
0.0 
0.0 
0.0 


RLT5 

745.30-^ 

866*416 

885^040 

903.681 

912.997 

931.630 

950.263 

959.579 

987.528 

1015.477 

1034.109 

1071.374 

1108.640 

1117.956 

L183.173 

1229.752 

1248.384 

1313.598 

1360. ISO 

1378.812 

1444.026 

1490.608 

1509.241 


FHi:3: 

a.o 

0.360 
0..361 
0;.362 
2.536 
2,542 
2.548 
L2.3  8  8 
12.427 
L2-.464 
6.244- 
6.266 
6.  .2  8  6 
L..480- 
L.A87 
1.4-92 
0.373 
0.375 
0.376 
0.376 
0.377 
0.3L7a 


XLO 

KLTl 

PHIl 

PLT2 

P}M2 

400.000 

0.0 

0.0 

0.0 

a.-o 

465.003 

0.0 

0.0 

0.0 

a..o 

475. COO 

0.0 

0.0 

Q.O 

o:..o: 

485.000 

0.0 

Q.O 

Q.o: 

a.,o 

490. J30 

0.0 

0.0 

0.0 

a..o 

500.000 

187.500 

1.082 

0.0 

a..o 

510. COO 

382.778 

1.082 

0.0 

a.o 

515.000 

433.688 

5.256 

0.0 

Q..0 

530.000 

544.064 

5.256 

0.0 

a..o 

545.000 

625.860 

5.256 

0.0 

0.0 

555.000 

672.289 

2.623 

0.0 

Q.O 

575.000 

753.476 

2.628 

0.0 

0.0 

595. 030 

324.555 

2.628 

0.0 

0.0 

600.000 

841.200 

0.618 

0.0 

0.0 

635.000 

948.366 

0.618 

0.0 

0.0 

660.000 

1013.887 

0.618 

0.0 

0.0 

670.000 

1045.730 

0.155 

0.0 

0.0 

705.000 

1135.654 

0.155 

0.0 

0.0 

730.033 

1196. 856 

0.155 

Q.O 

0.0 

740.000 

1220.771 

0.155 

O.Q 

0.0 

775.000 

1302.382 

0.155 

0.0 

0.0 

&00.00D 

1359.016 

0.155 

0.0 

0.0 

810.000 

1581.344 

0.0 

Q.O 

a..o 
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HEIGHT  OFF  THE  BOTTOM  = 


1197^., aOO  METERS: 


XLO 

RADIUS 

PHI 

RLT3 

PHI  3: 

400.003 

0.0 

0.0 

1098.289 

o:.o 

465.000 

0.0 

Q.O 

1276.761 

0.,196 

475.000 

0.0 

Q.O 

1304.219 

0..197 

485.0  00 

0.0 

a.o 

1331.676 

0.-197 

490.000 

0.0 

Q.O 

1345.404 

i:.3  84 

500.000 

0.0 

a.o 

1372.862 

1^339 

510.000 

0.0 

a.o 

1400.319 

U393 

515.000 

0.0 

Q.O 

1414..048 

6..7T5 

530.000 

0.0 

a.o 

L45  5.233 

6.,803_ 

545.000 

0.0 

Q.O 

L496..419 

6-.829 

555.000 

0.0 

0.0 

1523.876 

3:^.423 

575.000 

0.0 

0.0 

1578.791 

3.  .43  8 

595.000 

0.0 

a.o 

1633-705 

3i.453 

600.000 

0.0 

Q.O 

1647.434 

0.813 

635.000 

0.0 

Q.O 

1743.534 

0:..318 

660.000 

0.0 

Q.O 

18LZ.177 

0.822 

670.000 

0.0 

Q.O 

1839.635 

0.206 

705.000 

0.0 

0.0 

1935.735 

0.207 

730.000 

0.0 

O.Q 

2004.378 

0.207 

740.000 

0.0 

Q.O 

2031.835 

0...208 

775.000 

0.0 

Q.O 

2127.936 

0.203 

800.000 

0.0 

O.Q 

2196.579 

0.209 

810.000 

O.Q 

Q.Q 

2224  .-03.6 

OL..a 

XLO 

RLTl 

PHIL 

RLT2 

PHr2 

400.000 

0.0 

O.Q 

G..0 

Ou..O 

465.000 

0.0 

Q.Q 

0..0 

a.D 

475.000 

0.0 

Q.O 

0.0 

a.o: 

485.000 

0.0 

Q.O 

O.Q 

0.0 

490.000 

0.0 

0.0 

0.0 

O.D 

500.000 

0.0 

Q.O 

0.0 

0...0 

510.000 

527.833 

0.338 

0.0 

0..0 

515.000 

bll.460 

1.643 

0.0 

o:.o 

530.000 

734.567 

1.643 

0.0 

o.,o 

545.000 

909.296 

1.643 

0..0 

0...0 

555.000 

979. 561 

0.821 

0.0 

0.0 

575.000 

1101.499 

0.8  21 

0.0 

0.0 

595. 0  00 

1207.712 

0.821 

0.0 

o.o 

600.000 

1232.532 

0.193 

0.0 

OL.X) 

635.000 

1392.713 

0.193 

0.0 

0.0 

660. 000 

1496.637 

0.193 

O.Q 

0.0 

670.000 

1536.433 

0.048 

0.0 

0.0 

705.000 

1669.643 

0.048 

0.0 

0.0 

730.000 

1760.222 

0.048 

0.  0 

0.0 

740.000 

1795.598 

0.048 

0.0 

0.0 

775.000 

1916.269 

0.048 

O.Q 

0..  0 

800.000 

1999.965 

0.048 

O.Q 

O-.jO 

810.000 

2032.955 

Q.O 

Q..Q 

o:..o 
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COMPUTER  PKOGRAM 

IMPLICIT    REAL^a     (A-HtQ-Z) 

DIMEriSICN    Y(2,5)  .  THTA  (AJ  ,A  (4  ,47 ,.XLO  (  25  )  ,  ZPRT  (20)  t 
ISZ I Z5r 10) tSXLQ(25)  t  SY2 ( 25 , 1 0 )  ,SPH I( 25 , 1 0 )  , 

2SRLT1 (25»10) , S PHI L { 2? , 10 ) , S HLT2 ( 25 t IC) , SPH 1 2 ( 2 5 t 1 0 ) t  ■ 
4SRLT3{Z5,1J)  ,  SPH  13  (  25  t  U)  t  F  Y  (2  ,  5  )  ,  AZDEL  ( 1 } 

FltZZ,  YIrY2,DVT)  =  (  (  {  Z^GAfV  {  Y2'^  Y2  )  )^(Y2-{  3=^RTDEL/YI) 
l*tGA:VY2-ZZ)  )/  (3- {  0-ZZ)-»-2=-:=GAM/Y2  )  )+{  I  +  (  3  ^:=PTDHL/ Y  1 ) 
2=^tGAi1/{  Y2-Y2)  +  (GAM/-Y2-ZZ)-DVT/YI))  )/(3-(D-ZZ)  + 

32-ga:i/Y2) 

F2iZZrYirY2)=(YZ-Y^-(3*RTDEL-(GAM-(2LZ=;^y21  )/Yl)  )/{  3*Y2 

i-(D-ZZ)-i-2^GAM) 
C      THIS  PROGRAM  USES  THE  OPTIMAL  RUNGE"  KUTTA:  INTEGRATION 
C      METHOD. 

REAO{  5,  101)  GRAVtBETA  tALPH  ,  0  t-RTOEL  ,GAM  ,  XAQ- 
101  FORMAT Ca El  0.3) 

READ (5  ,101)  ZMAXtZO 
C      ZO  IS  THE  DEPTH  AT  WhFCH  THE  SCTURCE  IS  LOCATED.. 
C      ZMAX  IS  THE  SURFACE  DEPTH— HERE  IT  IS  ZERO^, 

READ  (5tlQl)  (THTA{ J)  ,J  =  1,4) 

READ (St  10  I)     ( (A( I, J)  , J  =  l,4) ,I=lt4) 
C  THE    'A*    AiJO    «THTa'     ARRAYS    ARE   USED    FOR    THEERUNGEi' 

C  KUTTA     INTEGRATION^ 

WRITEt6,666) 
666  FOR^'A~(»  1'  ) 

WRITEl  6,630)  &RAV,BETA,^LPH,a,^-TDEUGAK,AAC: 

630  FORMAT  <  ♦O'  ,7  (  5X  ,  EIQ  .3  ). ) 
WRITE (6,633)  ZMAXrZO 

633  FGR:-"aT(15X,'  ZMAX  =  •,Ei0.3,'  ZO  =  ",E10.3) 
WRITEi6T634)  {THTA{  J)  ,  J  =  lt4-) 

634  FGR;;AT{4(5X,E10.3)  ) 
KRITE(6,635)  {  {  A(  l,J)  ,  J  =  l,4)  »I.=  Lt4) 

635  F0Rr'1AT(4£15.3) 

READ  (.5, 130)  NZDELrNSKP,NXLOTNZPRT 
C      NZOEL  IS  THE  NUMBER  CF  INCREMENTS  FOR  INGEGR ATT  ON. . 
C      NSKP  IS  A  DEGUGGING  VARIABLE  USED  TO  LIMIT  THE 
C      LINES  OF  OUTPUT. 

C      riXLO  IS  THE  NUM6ER  OF  SUBBLES  QF  Dl FFEREIJT  INITIAL. 
C      RADII  TG  BE  ANALYZED.. 

C      MZPRT  IS  THE  NUMBER  C(^  DEPTHS.  FUR  WHICH- COMPARI  SON 
C      DATA  IS  TO  3E  GUTPUTTED. 
130  FQRMAT(3I10) 

READ (5, 101)  (AZDEL(J) ,J=1,NZDEL) 
C      AZDEL  IS  THE  INCREMENT  FOR  INTEGRATION.. 
WRITE  (6,  800)  NZDEL,.NSKP,NXLO,,NZPRT 
800  FCRMAT(15X,*  NZDEL  =  •,L6,«  NSKP  =  •  ,  1.8 1  •' '  NXLO  =- •  ,  ■ 
118,'  NZPRT  =  • ,13) 
WRITE(6r802)  (AZD£L(J) ,J=1, NZDEL) 
802  FORMAT (•  THE  AZDEL  APE  ',/,3E15.3) 
READ(5,101)  (XLO (J) ,J=1,NXL0) 
C      XLO  IS  THE  INITIAL  BUBBLE  RADIUS.. 
WRITE (6,631)  {XLG( J) rJ=ltNXLO) 

631  FORMAT  CO', 'THE  INITIAL  RADIUS  AR  E  •  ,/ ,  5X  ,  10E12.  3  ) 
READ  15, 101)  (ZPRT( J)  ,J  =  1,NZPRT) 

C  ZPRT     IS    THE    DEPTH    FOR    ViHICH    COMPARISON    DATA     IS    TO    BE 

C  OUTPUTTED. 

WRITE(6,632)  (ZPRT{ J) ,J=1, NZPRT) 

632  FORMAT  (»0','THE  INPUT  HEIGHTS  ARE ' , / » 5Xt IDE  12. 3) 
GB=GRAV/SETA 

C      NOW  COMPUTE  THE  INITIAL  BUBBLE  VELOCITY.. 

DC  80  K=I,NXL0 

DD2=DSQRT{XL0(K) ) 

DD3=RTDEL/DSuRT(GE/XAC) 

CALL    PHIS(XL0(K) ,PHIC) 

VT10  =  GB=^XL0(  K)  '"XL0(  K)  /ALPH 

VT20=DSQRT(G3-XL0{K)/XAC) 

VT0=  I  (-l.U+(  1.0+(4.0*VT10=:=VT10  )/  (  VT20?^VT20 )  ) 

i*=^0. 500000)/ 
2e(2,0-VT10^VT10)/{,VT2Q=^\/T20:)  )  )-VT10 

V0=VT0 
C  GB     IS    GRAV    DIVIDED    BY    BETA.. 

SXL0(K)=XL0{K)=^.1D7 
C      ALL  VARIABLES  BEGINNING  WITH  S  ARE  STORAGE  ARRAYS 
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C  FOR    THE    FINAL    OUTPUT. 

C  NOTE    THAT    RAOIL    ARE   READ    IN    AS  METERS    BU.r  OUTPUTTED 

C  IN    MICRONS. 

Y( 1, 1J=0. 

V(2,L)=XLQtK) 

IFLG=0 

DO    45     1=1,2 

DO    46    J=1t4 

FY<ItJ)=0. 
46    CGNTUiUE 
45    CCNTRiUE 

*>JRITEI6,636J     GB  ,  VTO  r  VO  »  VTIO  t  Vr20 

636  FGRMATl'l*  t2X, 'GB    =    ',£12.5,'     "^TO.      =    ••,£12-.  5:,  •• '  VO.  =-- •  ', 
lEi2.5, 

•    2*  VT13    =    »  ,E12*5r'  VT20    =    '  ,ElZ-.5:). 

WRITE(6»637)     ( Y ( J , 1) t J =L t 2 ) 

637  FORMAT{10X,613XtE12.5} J 

z=za 

KK=i 
KKK-0 

HHT=3.a 

C      HHT  15  THE  HEIGHT  OFF  THE  BOTTOM IT  I'S  taUAL.  TD  I-IO. 

ISO  =  0 
VT=VTO 
12=0 

ZTLDY=ZO 

KRITE(6,603)  VOr SXL3(K) tZO tZTLDYtVTO 
600  FORMATl '0»,/////T20Xr«  VG  =  ',£10. -3,'   METERS/SEC   ', 
120X, 

2»    LO    =    •,E13.3t'       microns    'tl,ZOXr*'    ZO    =^    ",FI0.3,- 
3«    METERS    «t/t 

420X,«     ZTLDY    =    'tFIQ.S,"    METERS.    ••,/',20X,  '     VO    TERM  INAL    = 
5    •, 

6Ei0.3,«    METERS/SEC ») 
WRITE{6,c66) 
11  =  0 
15=0 

XXL0=XL0lKi*.LD7 
701    11=11+1 

ZDEL=AZDELtlJ 
IS=IS+1 

IFdFLG.EQ.ll       GO    TO    661 
99    DC    23    J=l,4 

VTi=G3*Y(2,J)~Y(2T J)/ALPH 
VT2=DSQRT(G2=;=Y{2tJ) /XAC) 
VTR2=VTl*VTl/( VT2^VT2) 

VT=(  {-1.+DSQkT(1.  +  4.=^VTR2)  )  /  (  2 .  =!=VTR2)  )  ^!=VT1 
C      DVT  IS  DVT/DL  {A  DERIVATIVE). 

DVT=(3^-VTi/Y(2,  J)  )  /  (1  +  4'!^VTR2  )=^-0  .  500000-VT/Y{  2  ,  J  ) 
C      THE  INTEGRATION  IS  DONE  HERE.. 

FY(1tJ)=F1IZ+THTA( J)-ZDELt VT,Y(2tJ) tDVT) 
FY(2,J)=F2(Z+THTA(  J)-nZDEL,  VT,Y{2tJ  )) 
DO  2  5  N=1t2 

Y(N,  J  +  l)=Y{Nt  l)  +  {A{Jrl)^FY(N,l)  +  A(  JT2)=^FY{N,2)+Ai  J,3} 
I-FY(N,3J+ 
2AtJt4]-FY{Nr4) )*ZCEL 

rJ)  IS  LESS  THAN  ZERO,  IFLG  IS 
IfiTEGRATING  STEPS  ARE  SKIPPED 
PROGRAM. 
GO  TO  661 


c 

WHC;4  THE  RADIUS  Y(  2 

c 

SET  TO  ONE  AND  THE 

c 

FOR  THE  REST  OF  THE 
IF(Y{N,J+l).LT.O.) 

25 

CONTINUE 

23 

CONTINUE 
Y{1,1}=Y{1,5} 
Y(2,ii=Yt2,5) 
GO  TO  662 

661 

Yd, 11=1. 
Y{2,li=0. 
IFLG=1 
PHI=0. 

662 

Z=Z+ZDEL 

IF(Z.GT..iOOD-01) 

HHT=HHT+ZDEL 

GO  TO  999 
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IF( IS.EQ.MSKP)   GO  TO  28 

GO  TO  701 
28  IS=0 

IF(  IFLG.EQ.l)   GO  TO  663 

VT1  =  GB^Y{2t1)-Y(2t  D/ALPH 

VT2=DSQRT(GB-Y(2,1)/XAC) 

VTR2  =  VTl^-VTl/(  VT2='^VT2) 

VT  =  (  (-l.+DSgRT(l.  +  4.^^VTR2}  )  /  (  2  . '^VTR2  )  } -VTl 

GAM77  =  DEXP{Y(ia)  ) 

PHI  =  VTO'.--PHIO/(  VT^=GAM77) 
C      PHI  IS  THE  GENERAL  CASE  DISTRIBUTION. 
663  RZ=(D-Z)/(D-ZO) 

DD11=(D-Z0)/(D-Z) 

DDl=DSQRT( {D-ZO)/{D-Z) ) 

RLTl  =  (DDll*XLa(K>'^XL0(K)=i=XL0(K)  +  1.5*ALPH*RTDEL 
1*IZ-Z-Z0-ZG)/(GB^(0-Z) )) 

IF(RLT1.LT..1D-18)   RLT1=0. 

RLTl  =  RLTl--i^=i=  0.333  333333333 

PHI1=PHI0-RZ 
C      PHIL  IS  FOR  THE  VISCOUS  DRAG  MODEL. 

IF{RLT1.EQ.0.0)   PHI1=0.0 

RLT2  =  (  D01^'=DD2*DD2-DD2  +  D03=5'(  (D-Z}*{  1.-DD1*DD1*DD1 ) 
l*D=:^{l.-DDi)  )  ) 

IF(RLT2.LT.. lD-09)   FLT2=0^ 

RLT2  =  RLT2''=-0.666  666  6  66  666 

PHI2  =  PHI0---DSQRT(RZ) 
C      PHIZ  IS  FOR  THE  TURBULENT  DRAG,  SURFACE  ACTIVE 
C      MATERIALS  MODEL. 

IF(RLT2.EQ. 0.0)   PHI2=3.0 

RZ3  =  RZ*^^0.  333333333333 

RLT3=XL0(K)/RZ3 

VT13=GB-RLT3-RLT3/ALPH 

VT23=DSQRT<GB'^RLT3/XAC  ) 

VTR3  =  VT13-VT13/{VT23=!^VT23) 

VT3=(  (-l.+DSQRT{  1.+4.-VTR3)  )  /  (  2 .  =!--VTR3}  )'^VT13 

PHI3  =  PHI0^=VT0':'^RZ3/VT3 
C      PHI3  IS  FOR  THE  NO  GAS  DIFFUSION  CASE. 

Y2=Y(2,1)"1.D6 

XXXX=DABS(HHT-ZPRT(KK) ) 

RRLT1=RLT1*.1D7 

PRLT2=RLT2--'^.  ID7 

RRLT3=RLT3'-^.1D7 

IF(XXXX.GT.{0.1-ZDEL) )   GO  TO  223 

KKK=KKK+1 

SZ(K,KKK)=Z 

SY2(K,KKK)=Y(2,1)^-;'.1D7 

SPHI (K,KKK)=PHI 

SRLTl(KtKKK)=RLTl^^.lD7 

SPHI 1(K,KKK)=PHI I 

SRLT2(KtKKK)=RLT2-.1D7 

SPHI2(K,KKK)=PHI2 

SRLT3(  K,K.KK)  =RLT3-.1D7 

SPHI3(K,KKK)=PHI3 

KK=KK+1 
223  GO  TO  701 
999  CONTINUE 
80  CONTINUE 

DO  81  KKK=1,NZPRT 

WRIT5(6,82) 

82  FORMAT ( ' 1' ) 
WRITE(6,83)  ZPRT{KKK) 

83  FORMAT ( '0'  ,/////»30X, »  HEIGHT  OFF  THE  BOTTOM  =  •» 
1F13.3t '  METERS  ',///) 

WRITE(6,90) 
90  F0RMAT(23X,'        XLO   '»'       RADIUS't*        PHI 
I'        RLT3  • » •        PHI3  • ) 
DO  85  K=l,NXLO 

WRITE (6, 36)   SXLO{K) ,SY2(K,KKK) ,SPHI(K,KKK), 
1SRLT3( KtKKK)TSPHI3(K,KKK) 
85  CONTINUE 

V.'RITE(6,92) 
92  FORMAT(«0'  ,////) 
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WRIT 
91    FORM 
1» 

DO    8 

WRIT 

ISRLT 

86  FORM 

87  corjT 

CONT 
STOP 
END 


81 


E(6t91 
AT(23X 

RLT 
7  K  =  lt 
E(6,86 
2(KtKK 
AT(23X 
INUE 
INUE 


XLC   •,« 

PHI2« ) 


RLTl 


) 

2»,« 

NXLO 

)   SXLO(K),SRLTl(K,KKK),SPHIl(K,KKK)t 

K)  ,SPHI2(KtKKK) 

,5F11.3) 


PHIl  • 


9 

10 


SUBROUTINE  PHIS(R,PHTO) 
IMPLICIT  REAL=^8  (A-HtO-Z) 
R  IS  MEASURED  IN  METERS 


IF(R.LT. .46500-03) 

IF( (R.GE. .46500-03) 

IF( (R.GE. .49000-03) 

IF( (R.GE. .51500-03 ) 

IF{ (R.GE. .55500-03) 

IF( (R.GE. .60000-03) 

IF( (R.GE. .67  000- 03) 

IF( (R.GE. .74000-03) 

IF(R.GE. .81000-03) 

PHIO=0. 

GO  TO  10 

PHI0=1 . 

GO  TO  10 

PHI0=7. 

GO  TO  10 

PHI0=34. 

GO  TO  10 

PHI0=17. 

GO  TO  10 


GO  TO  1 

AMD. (R.LT. 
A^JO.  (R.LT. 
AND. (R.LT. 
AND. (R.LT. 
ANO. (R.LT. 
AND. (R.LT. 
AND. (R.LT. 
GO  TO  9 


.49000-03) ) 
.51500-03) ) 
.555000-03) } 
,60000-03) ) 
.67000-03) ) 
.74000-03) ) 
.8100D-03) ) 


GO  TO  2 
GO  TO  3 

GO  TO  4 
GO  TO  5 
GO  TO  6 
GO  TO  7 
GO  TO  8 


PH10=4 

GO  TO 

PHI0=1 

GO  TO 

PHI0=1 

GO  TO 

PHio---o. 

CONTINUE 

RETURN 

END 


10 
10 
10 
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